function [theta_NPD,fval_NPD,flag_NPD] = estimation_NPD_cvx(data_trans,J,N,m_regr,A_ineq,A_eq,b_ineq,b_eq,...
                                                            combos_all_unique,ntheta)

B_s = cell(J,1);
B_IV = cell(J,1);
w_mat = cell(J,1);

for j=1:J
   
    B_s{j} = single(zeros(N,size(combos_all_unique,1)));
    B_IV{j} = single(zeros(N,size(combos_all_unique,1)));
    w_mat{j} = single(zeros(size(combos_all_unique,1),size(combos_all_unique,1)));
        
end

for j=1:J
       
    B_s{j} = B_s_fun_4(data_trans{j},m_regr,combos_all_unique,N,J);        % (T,size(combos_all_unique,1)) array
    
    B_IV{j} = B_s{j};%B_IV_fun_3(data_trans{j},m_regr,combos_all_unique,N,J);
    
    temp = B_s{j}'*B_s{j};
    
    if det(B_s{j}*temp*B_s{j}')==0
        
        w_mat{j} = inv(B_s{j}'*B_s{j} + 0.01*eye(size(temp)));
        
    else w_mat{j} = inv(B_s{j}'*B_s{j});
    end
    
end

fun_criterion = @(x) estimation_criterion(x,data_trans,B_s,B_IV,w_mat,J);

ub = 100;
lb = -100;

if isempty(A_eq)==0
    
    cvx_begin quiet
    
    cvx_precision default
    variable theta_temp(ntheta,1)
    
    x = (reshape(theta_temp,[numel(theta_temp),1]));
    
    minimize( fun_criterion(x) )
    subject to
    A_ineq*x <= b_ineq
    A_eq*x == b_eq
    x<=ub
    x>=lb
    
    cvx_end
    
else    % i.e. if A_eq = []
    
    cvx_begin quiet
    
    cvx_precision default
    variable theta_temp(ntheta,1)
    
    x = (reshape(theta_temp,[numel(theta_temp),1]));
    
    minimize( fun_criterion(x) )
    subject to
    A_ineq*x <= b_ineq
    x<=ub
    x>=lb
    
    cvx_end
    
end

theta_NPD = theta_temp;
fval_NPD = cvx_optval;
flag_NPD = (strcmp(cvx_status,'Solved'));

end